function [n1, n2, n3] = demandModel(p1, p2, p3, N);
PI = 3.142;

m1 = 35;    %mean data price
m2 = 45;    %mean vedio price
m3 = 65;    %mean for bundle price

sd1 = 20;  %stand deviation of data willings to pay
sd2 = 30;  %stand deviation of vedio willings to pay
sd3 = 40;  %stand deviation of bundle willings to pay

r12 = -0.5;   %corelation number between data and vedio
r13 = -0.5;   %corelation number between data and bundle
r23 = -0.5;   %corelation number between vedio and bundle

w1 = 70;    %maximum willingnesss to pay for data
w2 = 80;    %maximum willingnesss to pay for vedio
w3 = 100;    %maximum willingnesss to pay for bundle

const1 = 2*PI*sd1*sd2*sd3*sqrt((1 - r12*r12)*(1 - r13*r13)*(1 - r23*r23));
const1 = 1/const1;

const2 = 2*(1 - r12*r12)*(1 - r13*r13)*(1 - r23*r23);
const2 = -1/const2;

n1 = 0;
n2 = 0;
n3 = 0;

% x = 1:1:n;
% y = 1:1:n;

for x1=1:1:w1
    for x2=1:1:w2
        for x3=1:w3;
            y1 = x1 - m1;
            y2 = x2 - m2;
            y3 = x3 - m3;

            part1 = y1*y1/(sd1*sd1);      	
            part2 = 2*r12*y1*y2/(sd1*sd2);
            part3 = y2*y2/(sd2*sd2);
            part4 = y3*y3/(sd3*sd3);
      	    part5 = 2*r13*y1*y3/(sd1*sd3);
            part6 = 2*r23*y2*y3/(sd2*sd3);
            
            P(x1,x2,x3) = const1*exp(const2*(part1-part2+part3+part4-part5-part6));   
            val = maxVal(x1-p1, x2-p2, x3-p3);
            if( val == 1)
                n1 = n1 + P(x1,x2,x3)*N;
            elseif (val == 2)
                n2 = n2 + P(x1,x2,x3)*N;
            else
                n3 = n3 + P(x1,x2,x3)*N;
            end
        end
    end
end
a = sum(sum(sum(P)))
% mesh(x,y,P)

function val = maxVal(v1, v2,v3);
if( v1 > v2)
    if( v1 > v3)
        val = 1;
    else
        val = 3;
    end
else
    if( v2 > v3 )
        val = 2;
    else
        val = 3;
    end
end
